R Notebook

if (R.version$major < 4) {
  warning("Not using R 4.0. It should.")
}
if (!require("pacman")) install.packages("pacman")

Loading required package: pacman
Warning: package ‘pacman’ was built under R version 4.3.2

pacman::p_load(
  tidyverse, # 
  conflicted,
  lubridate,
  simmer,
  simmer.plot,
  svglite
)

library(tidyverse)
library(conflicted)
library(lubridate)
conflicted::conflict_prefer("filter", "dplyr", c("base", "stats"))

[conflicted] Removing existing preference.[conflicted] Will prefer dplyr::filter over base::filter and stats::filter.

library(simmer)
library(simmer.plot)
conflicts_prefer(simmer::rollback)

[conflicted] Will prefer simmer::rollback over any other package.

conflicts_prefer(simmer::now)

[conflicted] Will prefer simmer::now over any other package.

conflicts_prefer(dplyr::select)

[conflicted] Will prefer dplyr::select over any other package.

options(
  repos=c(CRAN="https://cran.radicaldevelop.com/"),
  tidyverse.quiet = TRUE
)
# this will be substituted below and literally just is here because simmer doesn't know how to implement pipes and I hate it here

env <- simmer("hospital") 

Processo

Um paciente que está em lista de espera para ser operado:

1. Entra em contacto com o serviço de cirurgia, vários dias antes da sua operação,

• Através de uma chamada telefónica ou
• Através do envio de uma mensagem escrita numa app.

2. Os serviços administrativos atendem as chamadas ou lêem a mensagem via app e, mediante

• a informação presente no sistema dada pelo médico assistente e
• a informação fornecida pelo paciente, o administrativo decide se agenda
• uma consulta presencial antes da cirurgia, ou
• uma consulta telefónica, ou
• não há necessidade de qualquer consulta.

Contudo, este procedimento necessita de melhorias. O objectivo deste estudo é optimizar a gestão dos contactos (tempo total do processo e total de chamadas perdidas)

Situação atual

1. Apenas 2 administrativos no serviço, ambos começam às 8h e acabam às 18h (10h de trabalho) de cada dia útil.

2. O tempo entre chegadas tem exp(1/40) (pq media de 40/h)

3. Tipo de contacto:

• prob de telefone = 0.85
• prob de app = 0.15

4. Tipo de consulta:

• 0.3 não é agendado consulta,
• 0.6 consula por telefone,
• 0.1 presencial

5. Tempo de

• Conversa antes da decisao é norm(4,1^2)
• Tomada de decisao do adminstrativo é norm(1,0.25^2)
• Tempo de agendamento SE NECESSÁRIO é norm(1,0.25^2)

6. 5o telefone desliga apos 5 minutos de espera (e fica chamada perdida)

library(simmer)
library(simmer.bricks)

Warning: package ‘simmer.bricks’ was built under R version 4.3.2

library(simmer.plot)
conflicts_prefer(simmer::rollback)

[conflicted] Removing existing preference.[conflicted] Will prefer simmer::rollback over any other package.

conflicts_prefer(simmer::now)

[conflicted] Removing existing preference.[conflicted] Will prefer simmer::now over any other package.

conflicts_prefer(dplyr::select)

[conflicted] Removing existing preference.[conflicted] Will prefer dplyr::select over any other package.

# this will be substituted below and literally just is here because simmer doesn't know how to implement pipes and I hate it here
env <- simmer("hospital") 

Resources: Server: O hospital com 2 administradores Queue: Fila de espera com prioridade de clients em chamada Não há vantgem em escolher a app quando há chamada de espera Manager: vai atribuir o tipo de consulta (e Vai alterar o tipo de cliente se ele for rejeitado na chamada?) Gerador: vai gerir novos arrivals (exp(1/40)) arrival: cada paciente (que leva uma trajectory) trajectory: a receita de cada arrival

# atributes need to be numeric
tipo_de_consulta.PRESENCIAL <- 1
tipo_de_consulta.TELEFONICA <- 2
tipo_de_consulta.NAO_AGENDADA <- 3
tipos_de_consulta <- c(tipo_de_consulta.PRESENCIAL, tipo_de_consulta.TELEFONICA, tipo_de_consulta.NAO_AGENDADA)
tipos_de_consulta.probs <- c(0.1, 0.6, 0.3)

tipo_de_contacto.TELEFONE <- 1
tipo_de_contacto.APP <- 2
tipos_de_contacto <- c(tipo_de_contacto.TELEFONE, tipo_de_contacto.APP)
tipos_de_contacto.probs <- c(0.85, 0.15)

tempo_atendido <- function() {
  mean_time_triagem <- rnorm(1, mean = 4, sd = 1)
  mean_time_decision <- rnorm(1, mean = 1, sd = 0.25)
  time_to_timeout <- mean_time_triagem + mean_time_decision
  
  if (get_attribute(env, "tipo_de_consulta") != tipo_de_consulta.NAO_AGENDADA) {
    time_to_timeout <- time_to_timeout + rnorm(1, mean = 1, sd = 0.25)
  }
  
  return(time_to_timeout)
}

intime <- function() {
  now(env)%% 24*60 -> n
  return(8*60 < n & n < 18*60)
}

pessoa <- trajectory("pessoa") %>%
  # set_attribute("times_phone_rejected", 0) %>% 
  # set_attribute("phone_rejected", 0) %>%
  set_attribute("tipo_de_consulta", \() sample(tipos_de_consulta, 1, prob = tipos_de_consulta.probs)) %>%
  set_attribute("tipo_de_contacto", \() sample(tipos_de_contacto, 1, prob = tipos_de_contacto.probs), tag = "contacto") %>%
  set_prioritization(\() if (get_attribute(env, "tipo_de_contacto") == tipo_de_contacto.TELEFONE) c(1, -1, FALSE) else c(0 , -1, FALSE)) %>% # PRIORITY, DROP PRIORITY, RESTART
  # log_(\() paste("Vou pra fila agora, estou a usar", ifelse(get_attribute(env, "tipo_de_contacto") == tipo_de_contacto.TELEFONE, "o telefone", "a app"))) %>%
  renege_in(\() ifelse(get_attribute(env, "tipo_de_contacto") == tipo_de_contacto.TELEFONE, 5, Inf), out = 
              trajectory() %>% set_attribute("phone_rejected", 1)
              # O botas não gostou
              # trajectory("pessoa_try_again") %>% 
              # set_attribute("times_phone_rejected", \() get_attribute(env, "times_phone_rejected") + 1) %>% 
              # # leave if out of schedule
              # leave(\() {ifelse(intime(), 0, 1)}) %>% 
              # # log_(\() "Nao fui atendido, vou voltar pra fila") %>% 
              # rollback(target="contacto", times = 1) # não funciona, mudar pra numero?
            ) %>%
  seize("administrador", 1) %>%
  renege_abort() %>%
  # log_(\() paste("Before timeout")) %>% 
  set_attribute("atendido", 1) %>% 
  timeout(tempo_atendido) %>% 
  # log_(\() "After timeout") %>% 
  release("administrador", 1)
  # log_(\() "Leaving")

pessoa

trajectory: pessoa, 10 activities
{ Activity: SetAttribute | keys: [tipo_de_consulta], values: function(), global: 0, mod: N, init: 0 }
{ Activity: SetAttribute | [contacto] keys: [tipo_de_contacto], values: function(), global: 0, mod: N, init: 0 }
{ Activity: SetPrior     | values: function(), mod: N }
{ Activity: RenegeIn     | t: function(), keep_seized: 0 }
  Fork 1, stop,  trajectory: anonymous, 1 activities
  { Activity: SetAttribute | keys: [phone_rejected], values: [1], global: 0, mod: N, init: 0 }
{ Activity: Seize        | resource: administrador, amount: 1 }
{ Activity: RenegeAbort  |  }
{ Activity: SetAttribute | keys: [atendido], values: [1], global: 0, mod: N, init: 0 }
{ Activity: Timeout      | delay: function() }
{ Activity: Release      | resource: administrador, amount: 1 }

set.seed(1)
# 50 replicas
envs <- lapply(1:50, function(i) {
  env <<- simmer("hospital") 
  env %>% 
    add_generator("pessoa", pessoa, from_to(8*60, 18*60, \() rexp(1, 40/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_resource("administrador", schedule(c(8*60, 18*60), c(2, 0), period = 24*60)) %>%
    run(24*60*5)
})

envs %>% simmer.plot::get_mon_arrivals(ongoing = T) -> arrivals 
envs %>% simmer.plot::get_mon_attributes() -> attributes
envs %>% simmer.plot::get_mon_resources() -> resources

arrivals %>%
  filter(start_time != -1) %>% # n sei pq e q isto acontece
  arrange(start_time) %>% 
  select(-activity_time) %>% 
  mutate(
    day = ceiling(start_time/(24*60)),
    start_time_day = start_time %% (24*60) %>% seconds_to_period(),
    end_time_day = end_time %% (24*60) %>% seconds_to_period(),
    # replication = replication %>% as_factor(),
  ) -> arrivals.df1
arrivals.df1 %>% filter(replication == 1) %>% mutate(start_time_day = start_time_day %>% round(3), end_time_day = end_time_day %>% round(3))

# turn what's above into a function so it can be done for 2a and 2b
# get_arrivals <- function(mon_attributes) {
#   mon_attributes %>%
#     filter(start_time != -1) %>% # n sei pq e q isto acontece
#     arrange(start_time) %>% 
#     select(-activity_time) %>% 
#     mutate(
#       day = ceiling(start_time/(24*60)),
#       start_time_day = start_time %% (24*60) %>% seconds_to_period(),
#       end_time_day = end_time %% (24*60) %>% seconds_to_period(),
#       # replication = replication %>% as_factor(),
#     )
# }
# 
# pprint_arrivals <- function (arrivals.df1) {
#   arrivals.df1 %>% filter(replication == 1) %>% mutate(start_time_day = start_time_day %>% round(3), end_time_day = # end_time_day %>% round(3))
# }

# analise de eficiencia
# numero de chamadas perdidas
# add day to atributes from arrivals.df1
attributes %>% 
  left_join(arrivals.df1 %>% select(name, replication, day, start_time, end_time, end_time_day), by = c("name", "replication")) %>% 
  rename(atribute_time = time, pessoa_start_time = start_time, pessoa_end_time = end_time) %>% 
  select(atribute_time, name, key, value, replication, day, pessoa_start_time, pessoa_end_time) -> attributes.df1
attributes.df1 

# pessoas e os seus tipos de chamada
attributes.df1 %>% 
  filter(key == "tipo_de_contacto") %>% 
  select(name, replication, value) %>% 
  mutate(value = ifelse(value == tipo_de_contacto.TELEFONE, "telefone", "app")) %>%
  rename(contacto = value) -> contactos.df1

# percentagem de chamadas perdidas
attributes.df1 %>% 
  filter(key == "tipo_de_contacto", value == tipo_de_contacto.TELEFONE) %>% 
  left_join(attributes.df1 %>% filter(key == "phone_rejected") %>% mutate(gave_up = key) %>% select(name, replication, gave_up), join_by("name", "replication")) %>% 
  mutate(gave_up = if_else(gave_up %>% is.na, T, F), replication = replication %>% as_factor) %>% 
  select(name, replication, pessoa_end_time, gave_up) -> phone_calls
phone_calls

# phone_calls rate per replication
phone_calls %>% 
  group_by(replication) %>% 
  summarise(phone_calls = n(), gave_up = sum(gave_up)) %>% 
  mutate(phone_calls_per_hour = phone_calls/10, gave_up_rate = gave_up/phone_calls) %>% 
  arrange(desc(gave_up_rate)) -> phone_calls_rate
phone_calls_rate

phone_calls_rate %>% 
  ggplot(aes(x = reorder(replication, gave_up_rate), y = gave_up_rate)) +
  geom_col() + 
  ylim(0, 1) +
  theme(axis.text.x = element_text(angle = 70, hjust = 1)) +
  geom_hline(aes(yintercept = min(gave_up_rate)), linetype = "dashed", color = "red") +
  geom_hline(aes(yintercept = max(gave_up_rate)), linetype = "dashed", color = "blue") +
  xlab("Replication") + 
  ylab("Gave up rate")

phone_calls_rate %>% 
  summarise(
    min = min(phone_calls_per_hour), 
    mean = mean(phone_calls_per_hour), 
    max = max(phone_calls_per_hour)
    )

# evolucao ao longo do dia DO QUE?
attributes.df1 %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

# get queue size histogram
resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,250)

NA
NA

# join these 2 plots
attributes.df1 %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,875)

  
attributes.df1 %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

attributes.df1 %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +

  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +

  theme(legend.position = "none")

# tempo de espera daqueles que concluiram, speardao por chamada e mensagem
attributes %>% 
  filter(key == "atendido") %>% 
  left_join(arrivals.df1 %>% select(start_time, name, replication), join_by("name", "replication")) %>% 
  left_join(contactos.df1, join_by("name", "replication")) %>% 
  mutate(time_waiting = time - start_time) %>% 
  select(name, value, replication, time_waiting, contacto) -> waiting_times
waiting_times

# histograma (n tou a perceber pq é q há waiting times maiores que 5)
waiting_times %>% 
  filter(contacto == "telefone") %>%
  ggplot(aes(x = time_waiting)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")

waiting_times %>% 
  filter(contacto == "app") %>%
  ggplot(aes(x = time_waiting/60)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")

# regular waiting time stats without graphs
waiting_times %>% 
  group_by(contacto) %>% 
  summarise(mean = mean(time_waiting), sd = sd(time_waiting), median = median(time_waiting)) 

waiting_times %>% 
  group_by(contacto) %>% 
  summarise(quantile_10 = quantile(time_waiting, 0.1), quantile_20 = quantile(time_waiting, 0.2), quantile_30 = quantile(time_waiting, 0.3), quantile_40 = quantile(time_waiting, 0.4), quantile_50 = quantile(time_waiting, 0.5), quantile_60 = quantile(time_waiting, 0.6), quantile_70 = quantile(time_waiting, 0.7), quantile_80 = quantile(time_waiting, 0.8), quantile_90 = quantile(time_waiting, 0.9), quantile_95 = quantile(time_waiting, 0.95), quantile_99 = quantile(time_waiting, 0.99))

resources %>% plot(metric = "utilization")

resources %>% plot(metric = "usage", steps = T)

Pergunta 2

Alíneas com mudanças relativas às interrupções

# mudanças
envs <- lapply(1:50, function(i) {
  env <<- simmer("hospital") 
  env %>% 
    add_generator("pessoa_de_manha", pessoa, from_to(8*60, 10*60, \() rexp(1, (120/2)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_generator("pessoa_de_tarde", pessoa, from_to(10*60, 16*60, \() rexp(1, (240/6)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_generator("pessoa_de_noite", pessoa, from_to(16*60, 18*60, \() rexp(1, (40/2)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_resource("administrador", schedule(
      c(8 , 9 , 9.5 , 12 , 15.5 ,18)*60, 
      c(  1 , 2,    4,   3,     1,  0), 
      period = 24*60
    )) %>%
    run(24*60*5)
})

envs %>% simmer.plot::get_mon_arrivals(ongoing = T) -> arrivals
envs %>% simmer.plot::get_mon_attributes() -> attributes
envs %>% simmer.plot::get_mon_resources() -> resources

2.b)

# b)
arrivals %>%
  filter(start_time != -1) %>% # n sei pq e q isto acontece
  arrange(start_time) %>% 
  select(-activity_time) %>% 
  mutate(
    day = ceiling(start_time/(24*60)),
    start_time_day = start_time %% (24*60) %>% seconds_to_period(),
    end_time_day = end_time %% (24*60) %>% seconds_to_period(),
    replication = replication %>% as_factor(),
  ) -> arrivals.df2b
arrivals.df2b %>% filter(replication == 1) %>% mutate(start_time_day = start_time_day %>% round(3), end_time_day = end_time_day %>% round(3))

arrivals.df2b$replication <- as.integer(as.character(arrivals.df2b$replication))

# tempo de espera daqueles que concluiram, speardao por chamada e mensagem
arrivals.df2b %>% 
  filter(finished == TRUE) %>% 
  mutate(waiting_time = end_time - start_time) %>% 
  left_join(attributes %>% filter(key == "tipo_de_contacto"), by = c("name", "replication")) %>% 
  mutate(tipo_de_contacto = ifelse(value == tipo_de_contacto.TELEFONE, "telefone", "mensagem") %>% as_factor, replication = replication %>% as_factor) %>%
  select(name, replication, waiting_time, tipo_de_contacto) -> waiting_times.2b
waiting_times.2b

# analise de eficiencia
# numero de chamadas perdidas
# add day to atributes from arrivals.df2b
attributes %>% 
  left_join(arrivals.df2b %>% select(name, replication, day, start_time, end_time, end_time_day), by = c("name", "replication")) %>% 
  rename(atribute_time = time, pessoa_start_time = start_time, pessoa_end_time = end_time) %>% 
  select(atribute_time, name, key, value, replication, day, pessoa_start_time, pessoa_end_time) -> attributes.df2b
attributes.df2b 

# pessoas e os seus tipos de chamada
attributes.df2b %>% 
  filter(key == "tipo_de_contacto") %>% 
  select(name, replication, value) %>% 
  mutate(value = ifelse(value == tipo_de_contacto.TELEFONE, "telefone", "app")) %>%
  rename(contacto = value) -> contactos.df2b

# percentagem de chamadas perdidas
attributes.df2b %>% 
  filter(key == "tipo_de_contacto", value == tipo_de_contacto.TELEFONE) %>% 
  left_join(attributes.df2b %>% filter(key == "phone_rejected") %>% mutate(gave_up = key) %>% select(name, replication, gave_up), join_by("name", "replication")) %>% 
  mutate(gave_up = if_else(gave_up %>% is.na, T, F), replication = replication %>% as_factor) %>% 
  select(name, replication, pessoa_end_time, gave_up) -> phone_calls
phone_calls

# phone_calls rate per replication
phone_calls %>% 
  group_by(replication) %>% 
  summarise(phone_calls = n(), gave_up = sum(gave_up)) %>% 
  mutate(phone_calls_per_hour = phone_calls/10, gave_up_rate = gave_up/phone_calls) %>% 
  arrange(desc(gave_up_rate)) -> phone_calls_rate
phone_calls_rate

phone_calls_rate %>% 
  ggplot(aes(x = reorder(replication, gave_up_rate), y = gave_up_rate)) +
  geom_col() + 
  ylim(0, 1) +
  theme(axis.text.x = element_text(angle = 70, hjust = 1)) +
  geom_hline(aes(yintercept = min(gave_up_rate)), linetype = "dashed", color = "red") +
  geom_hline(aes(yintercept = max(gave_up_rate)), linetype = "dashed", color = "blue") +
  xlab("Replication") + 
  ylab("Gave up rate")

phone_calls_rate %>% 
  summarise(
    min = min(phone_calls_per_hour), 
    mean = mean(phone_calls_per_hour), 
    max = max(phone_calls_per_hour)
    )

# evolucao ao longo do dia DO QUE?
attributes.df2b %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

# get queue size histogram
resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,250)

NA

# join these 2 plots
attributes.df2b %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,875)

  
attributes.df2b %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ggtitle("Tamanho da fila e chamadas rejeitadas 2b) por dia") +
  xlab("Tempo de simulação (horas)") + 
  ylab("Número de chamadas rejeitadas")

attributes.df2b %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +

  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  theme(legend.position = "none") +
  ggtitle("Tamanho da fila e chamadas rejeitadas 2b) nos 5 dias") +
  xlab("Tempo de simulação (horas)") + 
  ylab("Número de chamadas rejeitadas")

# tempo de espera daqueles que concluiram, speardao por chamada e mensagem
attributes %>% 
  filter(key == "atendido") %>% 
  left_join(arrivals.df2b %>% select(start_time, name, replication), join_by("name", "replication")) %>% 
  left_join(contactos.df2b, join_by("name", "replication")) %>% 
  mutate(time_waiting = time - start_time) %>% 
  select(name, value, replication, time_waiting, contacto) -> waiting_times
waiting_times

# histograma (n tou a perceber pq é q há waiting times maiores que 5)
waiting_times %>% 
  filter(contacto == "telefone") %>%
  ggplot(aes(x = time_waiting)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")+
  ylim(0, 1000)

waiting_times %>% 
  filter(contacto == "app") %>%
  ggplot(aes(x = time_waiting/60)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")

# regular waiting time stats without graphs
waiting_times %>% 
  group_by(contacto) %>% 
  summarise(mean = mean(time_waiting), sd = sd(time_waiting), median = median(time_waiting)) 

waiting_times %>% 
  group_by(contacto) %>% 
  summarise(quantile_10 = quantile(time_waiting, 0.1), quantile_20 = quantile(time_waiting, 0.2), quantile_30 = quantile(time_waiting, 0.3), quantile_40 = quantile(time_waiting, 0.4), quantile_50 = quantile(time_waiting, 0.5), quantile_60 = quantile(time_waiting, 0.6), quantile_70 = quantile(time_waiting, 0.7), quantile_80 = quantile(time_waiting, 0.8), quantile_90 = quantile(time_waiting, 0.9), quantile_95 = quantile(time_waiting, 0.95), quantile_99 = quantile(time_waiting, 0.99))

resources %>% plot(metric = "utilization")

resources %>% plot(metric = "usage", steps = T)

2.a) - - - - - - - - - - - -

# a) set_prioritization(\() if (get_attribute(env, "tipo_de_consulta") == tipo_de_consulta.TELEFONICA) c(1, -1, FALSE) else c(0 , -1, FALSE))
pessoa__ <- join(
  pessoa[1:2], 
  trajectory("change_prio") %>% 
     set_prioritization(\() if (get_attribute(env, "tipo_de_contacto") == tipo_de_contacto.TELEFONE) c(1, 1, FALSE) else c(0 , 0, FALSE)),
  pessoa[-(1:3)]
)
pessoa__

trajectory: pessoa, 10 activities
{ Activity: SetAttribute | keys: [tipo_de_consulta], values: function(), global: 0, mod: N, init: 0 }
{ Activity: SetAttribute | [contacto] keys: [tipo_de_contacto], values: function(), global: 0, mod: N, init: 0 }
{ Activity: SetPrior     | values: function(), mod: N }
{ Activity: RenegeIn     | t: function(), keep_seized: 0 }
  Fork 1, stop,  trajectory: anonymous, 1 activities
  { Activity: SetAttribute | keys: [phone_rejected], values: [1], global: 0, mod: N, init: 0 }
{ Activity: Seize        | resource: administrador, amount: 1 }
{ Activity: RenegeAbort  |  }
{ Activity: SetAttribute | keys: [atendido], values: [1], global: 0, mod: N, init: 0 }
{ Activity: Timeout      | delay: function() }
{ Activity: Release      | resource: administrador, amount: 1 }

envs <- lapply(1:50, function(i) {
  env <<- simmer("hospital") 
  env %>% 
    add_generator("pessoa_de_manha", pessoa__, from_to(8*60, 10*60, \() rexp(1, (120/2)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_generator("pessoa_de_tarde", pessoa__, from_to(10*60, 16*60, \() rexp(1, (240/6)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_generator("pessoa_de_noite", pessoa__, from_to(16*60, 18*60, \() rexp(1, (40/2)/60), every = 24*60, arrive = F), mon = 2) %>% 
    add_resource("administrador", schedule(
      c(8 , 9 , 9.5 , 12 , 15.5 ,18)*60, 
      c(  1 , 2,    4,   3,     1,  0), 
      period = 24*60
    )) %>%
    run(24*60*5)
})

envs %>% simmer.plot::get_mon_arrivals(ongoing = T) -> arrivals
envs %>% simmer.plot::get_mon_attributes() -> attributes
envs %>% simmer.plot::get_mon_resources() -> resources

arrivals %>%
  filter(start_time != -1) %>% # n sei pq e q isto acontece
  arrange(start_time) %>% 
  select(-activity_time) %>% 
  mutate(
    day = ceiling(start_time/(24*60)),
    start_time_day = start_time %% (24*60) %>% seconds_to_period(),
    end_time_day = end_time %% (24*60) %>% seconds_to_period(),
    replication = replication %>% as_factor(),
  ) -> arrivals.df2a
arrivals.df2a %>% filter(replication == 1) %>% mutate(start_time_day = start_time_day %>% round(3), end_time_day = end_time_day %>% round(3)) 

arrivals.df2a$replication <- as.integer(as.character(arrivals.df2a$replication))

# analise de eficiencia
attributes %>% 
  left_join(arrivals.df2a, by = c("name", "replication") ) %>% 
  select(name, key, value, replication, day, end_time, end_time_day) -> atributes.df2a
atributes.df2a %>% 
  filter(key == "phone_rejected")

# analise de eficiencia
# numero de chamadas perdidas
# add day to atributes from arrivals.df2a
attributes %>% 
  left_join(arrivals.df2a %>% select(name, replication, day, start_time, end_time, end_time_day), by = c("name", "replication")) %>% 
  rename(atribute_time = time, pessoa_start_time = start_time, pessoa_end_time = end_time) %>% 
  select(atribute_time, name, key, value, replication, day, pessoa_start_time, pessoa_end_time) -> attributes.df2a
attributes.df2a 

# pessoas e os seus tipos de chamada
attributes.df2a %>% 
  filter(key == "tipo_de_contacto") %>% 
  select(name, replication, value) %>% 
  mutate(value = ifelse(value == tipo_de_contacto.TELEFONE, "telefone", "app")) %>%
  rename(contacto = value) -> contactos.df2a

# percentagem de chamadas perdidas
attributes.df2a %>% 
  filter(key == "tipo_de_contacto", value == tipo_de_contacto.TELEFONE) %>% 
  left_join(attributes.df2a %>% filter(key == "phone_rejected") %>% mutate(gave_up = key) %>% select(name, replication, gave_up), join_by("name", "replication")) %>% 
  mutate(gave_up = if_else(gave_up %>% is.na, T, F), replication = replication %>% as_factor) %>% 
  select(name, replication, pessoa_end_time, gave_up) -> phone_calls
phone_calls

# phone_calls rate per replication
phone_calls %>% 
  group_by(replication) %>% 
  summarise(phone_calls = n(), gave_up = sum(gave_up)) %>% 
  mutate(phone_calls_per_hour = phone_calls/10, gave_up_rate = gave_up/phone_calls) %>% 
  arrange(desc(gave_up_rate)) -> phone_calls_rate
phone_calls_rate

phone_calls_rate %>% 
  ggplot(aes(x = reorder(replication, gave_up_rate), y = gave_up_rate)) +
  geom_col() + 
  ylim(0, 1) +
  theme(axis.text.x = element_text(angle = 70, hjust = 1)) +
  geom_hline(aes(yintercept = min(gave_up_rate)), linetype = "dashed", color = "red") +
  geom_hline(aes(yintercept = max(gave_up_rate)), linetype = "dashed", color = "blue") +
  xlab("Replication") + 
  ylab("Gave up rate")

phone_calls_rate %>% 
  summarise(
    min = min(phone_calls_per_hour), 
    mean = mean(phone_calls_per_hour), 
    max = max(phone_calls_per_hour)
    )

# evolucao ao longo do dia DO QUE?
attributes.df2a %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

# get queue size histogram
resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none")

resources %>% 
  mutate(day = ceiling(time/(24*60))) %>% 
  group_by(replication) %>%
  reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
  mutate(replication = replication %>% as_factor()) %>% 
  ggplot(aes(x = time, y = queue_size)) +
  geom_step(aes(color = replication), alpha = 0.3) + 
  geom_smooth(method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,250)

NA

# join these 2 plots
attributes.df2a %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ylim(0,875)

  
attributes.df2a %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  facet_wrap(vars(day), scale = "free") +
  theme(legend.position = "none") +
  ggtitle("Tamanho da fila e chamadas rejeitadas 2a) por dia") +
  xlab("Tempo de simulação (horas)") + 
  ylab("Número de chamadas rejeitadas")

attributes.df2a %>% 
  filter(key == "phone_rejected") %>% 
  mutate(day = day %>% as_factor(), replication = replication %>% as_factor) %>% 
  group_by(replication) %>% 
  reframe(time_rejected = pessoa_end_time/60, sumcum_per_day = cumsum(value), day = day) %>% 
  ggplot(aes(x = time_rejected, y = sumcum_per_day)) +
  geom_line(aes(color = replication), alpha = 0.3) +
  geom_smooth() +

  theme(legend.position = "none") +
  geom_step(data = resources %>% 
              mutate(day = ceiling(time/(24*60))) %>% 
              group_by(replication) %>%
              reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
              mutate(replication = replication %>% as_factor()), 
            aes(x = time, y = queue_size), alpha = 0.3) + 
  geom_smooth(data = resources %>% 
                mutate(day = ceiling(time/(24*60))) %>% 
                group_by(replication) %>%
                reframe(day = day %>% as_factor(), time = time/60, queue_size = queue) %>% 
                mutate(replication = replication %>% as_factor()), 
              aes(x = time, y = queue_size), method="lm") +
  theme(legend.position = "none") +
  ggtitle("Tamanho da fila e chamadas rejeitadas 2a) nos 5 dias") +
  xlab("Tempo de simulação (horas)") + 
  ylab("Número de chamadas rejeitadas")

# tempo de espera daqueles que concluiram, speardao por chamada e mensagem
attributes %>% 
  filter(key == "atendido") %>% 
  left_join(arrivals.df2a %>% select(start_time, name, replication), join_by("name", "replication")) %>% 
  left_join(contactos.df2a, join_by("name", "replication")) %>% 
  mutate(time_waiting = time - start_time) %>% 
  select(name, value, replication, time_waiting, contacto) -> waiting_times
waiting_times

# histograma (n tou a perceber pq é q há waiting times maiores que 5)
waiting_times %>% 
  filter(contacto == "telefone") %>%
  ggplot(aes(x = time_waiting)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")+ 
  ylim(0, 1000)

waiting_times %>% 
  filter(contacto == "app") %>%
  ggplot(aes(x = time_waiting/60)) +
  geom_histogram(bins = 100) +
  theme(legend.position = "none")

# regular waiting time stats without graphs
waiting_times %>% 
  group_by(contacto) %>% 
  summarise(mean = mean(time_waiting), sd = sd(time_waiting), median = median(time_waiting)) 

waiting_times %>% 
  group_by(contacto) %>% 
  summarise(quantile_10 = quantile(time_waiting, 0.1), quantile_20 = quantile(time_waiting, 0.2), quantile_30 = quantile(time_waiting, 0.3), quantile_40 = quantile(time_waiting, 0.4), quantile_50 = quantile(time_waiting, 0.5), quantile_60 = quantile(time_waiting, 0.6), quantile_70 = quantile(time_waiting, 0.7), quantile_80 = quantile(time_waiting, 0.8), quantile_90 = quantile(time_waiting, 0.9), quantile_95 = quantile(time_waiting, 0.95), quantile_99 = quantile(time_waiting, 0.99))

resources %>% plot(metric = "utilization")

resources %>% plot(metric = "usage", steps = T)